//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Declaracoes das funcoes da biblioteca c++ para programa de modelagem de
// curva de luz. Versao C++ da biblioteca roche.h para modelar curvas de luz
// da estrela secundaria em Variaveis Cataclismicas. Esta versao deve ser mais
// generica e permitir outras geometrias de forma mais simples que a versao
// C. Por exeplo, calcular curva de luz para diferentes geometrias, disco,
// estrela esferica, estrela preenchendo lobo de roche e estrela nao 
// preenchendo o lobo de roche. 
//
// Comentarios completos, com log de atividades no arquivo light_curve.cpp
//
// Tiago Ribeiro -- LinuxC++ -- UFSC -- 2006
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Log de atualizacoes:
//
// 18 Out 2006 - Primeira versao funcional, com o basico. Resultados 
//             equivalentes aos do programa em C.
// 30 Out 2006 - Primeira versao de SPOT_MAKER (nao testada ainda) e 
//             inclusao de mais Leis de Obscurecimento de borda; 
//             linear, quadratica, logaritmica e nao linear de Claret 2000.
// -- Mai 2007 - Implementando rotinas de pixelizacao para facilitar a criacao
//             dos mapas de distribuicao de brilho da secundaria com OpenGl. 
//             Pensando em implementar uma nova rotina de pixelizacao para o 
//             programa (Tegmark M. 1996).
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#ifndef LIGHT_CURVE_DEF

#define LIGHT_CURVE_DEF

#include <cstdlib>
#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <vector>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include "lib_zabot/aborta.h"
#include "histograma.h"
//#include "histograma.cpp"
#include "namespace_fis_const.h"

using namespace std;

struct PixMap {
  vector<int> npix; //numero de pixels por coluna
  vector<double> theta;
  vector<double> phi;
  vector<double> px; //posicao x
  vector<double> py; //posicao y
  vector<double> pz; //posicao z
};

struct PixOpenGL {
  double theta[3];
  double phi[3];
  double px[3];
  double py[3];
  double pz[3];
  double I;
};

struct surface {
  vector<double> theta;
  vector<double> phi;
  vector<double> px; //posicao x
  vector<double> py; //posicao y
  vector<double> pz; //posicao z
  vector<double> nx; //versor x
  vector<double> ny; //versor y
  vector<double> nz; //versor z
  vector<double> I;  //intensidade 
  vector<double> vx; //velocide x
  vector<double> vy; //velocide y
  vector<double> vz; //velocide z
};

struct pixel {
  double theta;
  double phi;
  double px; //posicao x
  double py; //posicao y
  double pz; //posicao z
  double nx; //versor x
  double ny; //versor y
  double nz; //versor z
  double I;  //intensidade 
  double vx; //velocide x
  double vy; //velocide y
  double vz; //velocide z
};


struct par_list {
  double i; //inclinacao
  double q; //razao de massa
  double f_fac; //filling factor (fator multiplicativo a potencial do L1)
  double obs_grav; //obscurecimento gravitacional
  double g_pole; //gravidade no polo
};


// Classe padrao para gerar curva de luz! 
// possui apenas as coisas basicas!

class Light_Curve {

  // definicao do vetor ->  vector<vector<double> > par; (definido mais abaixo)
  // um vetor pra cada tipo de parametro! beeem melhor! Use como quiser mas os 
  // 3 primeiros sao utilizados como default para os parametros orbitais, 
  // atmosfericos e obscurecimento, os seguintes sao utilizados pra definir os
  // os spots. O formato considerado eh como segue:
  // par[0][i]:
  // 0 - i
  // 1 - M_1
  // 2 - M_2
  // 3 - P_orb
  // par[1][i]:
  // 0 - fluxo da estrela em fase 0		       
  // 1 - coef. para obscurecimento gravitacional
  // 2 - lei de obscurecimento de borda (veja help na funcao!)  
  // 3 - numero de spots na superficie da estrela      
  // 4 - Temperatura no polo da estrela
  // 5 - Filling factor (Fator multiplicativo ao potencial no L1)
  // 6 - Albedo
// par[2][i]:
  // 0 a n - coeficientes de obscurecimento de borda
  // par[3][i] ate par[3+par[1][3]][i]:
  // 0 - posicao em theta
  // 1 - posicao em phi
  // 2 - dtheta
  // 3 - dphi
  // 4 - Is

  // Leis de obscurecimento de borda.
  double sqrt_law(const double nv);
  double lin_law(const double nv); // Lei Linear
  double quad_law(const double nv); // Lei quadratica
  double ln_law(const double nv); // Lei logaritmica
  double nlin_claret_law(const double nv);//Lei proposta por Claret A. 2000AA.363..1081
  bool _vorb;
  double _psize; // Guarda o tamanho do pixel
 protected:

  surface star;
  vector<vector<double> > par; // vetor para guardar os parametros do sistema!
  double default_ec(int,const vector<double> &);
  //double (*eclipse)(...);
  double eclipse(int npix, double (*ec_func)());
 public:
	Light_Curve(){_vorb=true;};//Inicializando valores internos
  Light_Curve(const vector<vector<double> > &data_par);
  virtual double flux(const double fase, bool eclipsa = false);
  int close_pix(const double, const double, const double, const double, double*& );
	//double flux(const double fase, const double *(*)(const vector<double> &,const vector<double> &));
  //double flux(const double fase, double (*ec_fnc)(void * PAR),void * ec_par);
  //template <typename T> double flux(const double fase, double (&ec_fnc)(...),T &ec_par);

  template <typename T> double flux(const double fase,T &ec_par, double (&ec_fnc)(vector<double>&,vector<double>&,T&));
	
  double obscurecimento_de_borda(const double nv, const int law);
  int npar_obs_borda(const int law);
  int size(void);
  void sum_flux_pix(const int n,const double f); //soma f ao pixel n
  void set_flux_pix(const int n,const double f);
	double pix_size(void){
		return _psize;
	};
	void pix_size(double val){
		_psize = val;
	};
	
  int proc_pix(const double theta, const double phi); //devolve posicao de 
                                                      //pixel mais 
                                                      //proximo de theta e phi.

  int proc_pix(const double x, const double y, const double z);
  void spot_maker(void);
  void spot_maker_test(void);
	void albedo(pixel&,int);
  surface Star(void);
  double PixVel(int pix, double fase); //devolve velocidade do pixel pix
  //devolve a posicao no histograma de formato format de pixel com 
  //velocidade vel
  int PixVelPos(double vel, const vector<double> &format);

  double VFunc(int n, double fase);

  Histograma rad_prof(const double fase, const vector<double> &format);

  Histograma rad_prof(const double, const vector<double> &, const vector<int> &);

  void inline pixmap(const double fase, Histograma &input, vector<int> &output)
    {
#ifdef VERBOSE
      if(output.size() != 0)
	cout << "# ATENCAO: Vetor de saida inicializado. Resetando vetor e continuando...\n";
#endif
      
      output.resize(0);

      //cout << "# [pixmap] -" << star.theta.size() << endl;

      for(int i = 0; i < star.theta.size(); i++)
	output.push_back( ( (VFunc(i,fase) > 0.0 )? input.PixPos(PixVel(i,fase)) : -1 ) );
      
    };
    void inline vorb(bool val){_vorb=val;};
    bool inline vorb(void){return _vorb;};
};

// Estrela preenchendo o lobo de roche!

class Roche: public Light_Curve {
  static const int N_PAR_ORBITAL = 4; // numero de parametros orbitais
  // descricao:
  // 0 - i
  // 1 - M_1
  // 2 - M_2
  // 3 - P_orb
  static const int N_PAR_ATMOSF = 5; // numero de parametros atmosfericos
  // descricao:
  // 0 - brilho da estrela em fase 0		       
  // 1 - coef. para obscurecimento gravitacional
  // 2 - numero de spots na superficie da estrela      
  // 3 - numero de coef. para obscurecimento de borda  
  // 4 - Temperatura no polo da estrela                

  struct par_list roche_par;
  double rl1; // raio do ponto lagrangiano interno
  double F20;
  int __Dec__;
  int _size_;
  
 public:
  Roche();
  Roche(const int dec, const vector<vector<double> > &data_par);
  Roche(char *pixfile, const vector<vector<double> > &data_par);
  Roche(const surface &StarSurf, const vector<vector<double> > &data_par);
  template <typename T> void NewI(T &);
  void NewI(const gsl_vector &newi);
  void NewI(const vector<double> &newi)
  {
    if(newi.size() != star.theta.size()) aborta("Nova mapa de tamanho diferente do esperado!");
    double  qq = 1./roche_par.q; // lobo da secundaria
    for(int i = 0; i < newi.size(); i++){
      double rlob = ache_r(star.theta[i], star.phi[i]);
      double x = star.px[i];// *rlob;
      double y = star.py[i];// *rlob;
      double z = star.pz[i];// *rlob;
      double nr = sqrt(x*x+y*y+z*z);
      double n_r = sqrt((1-x)*(1-x)+(y*y)+(z*z));
      nr = nr*nr*nr;
      n_r = n_r*n_r*n_r;
      double nx = -x/nr+(qq*(1-x)/n_r)+(qq+1)*(x)-qq;
      double ny = -y*(1./nr+qq/n_r-(qq+1));
      double nz = -z*(1./nr+qq/n_r);
      double np = sqrt(nx*nx+ny*ny+nz*nz);
      //star.I[w] = pow(np/roche_par.g_pole, roche_par.obs_grav);

      star.I[i]=newi[i]*pow(np/roche_par.g_pole, roche_par.obs_grav);
    }
  };
  void CalcPixSize(void);
  void openGL_Image(const int dec, surface &Image);
  void openGL_Image_polar(const int dec, surface &Image);
  void GL_pix(const int dec, surface &Image);
  void GL_pix_polar(const int dec, surface &Image);
  double flux2(void) {return F20;};
  void flux2(const double setf2){F20 = setf2;};
  void set_RL1(const double setrl1){rl1=setrl1;};
  double inline get_RL1(void){return rl1;};
  par_list &rpar(void){return roche_par;};
  void rpar(const par_list &set_rp)
  {
    roche_par.i=set_rp.i;
    roche_par.q=set_rp.q;
    roche_par.f_fac=set_rp.f_fac;
    roche_par.obs_grav=set_rp.obs_grav;
    roche_par.g_pole=set_rp.g_pole;
  };
  double Phi0(void) {
    return 2./((1+roche_par.q)*rl1)+2.*roche_par.q/((1+roche_par.q)*sqrt(pow(rl1,2.)-(2.*rl1)+1))+pow((rl1)-(roche_par.q/(1+roche_par.q)),2.);
  };
  double inline Phi(double xx,double yy , double zz)
  { 
    double &q = roche_par.q;
    double R = sqrt(xx*xx+yy*yy+zz*zz);
    double r_2 = sqrt(pow(R,2.)-2.*xx+1);
    double a = 2./((1+q)*R);
    double b = 2.*q/((1+q)*r_2);
    double c = pow(xx-(q/(1+q)),2.);
    double d = pow(yy,2.);
    return a+b+c+d;
  };
  int inline npix(void){return star.theta.size();};
  void healpixmap_2_vtk(int dec, ostringstream &vtk_out);
  double gpole(void);
  void lagrangeana(void);
  double ache_r(double v0, double u0);
  double ache_r2(double v0, double u0);
  void v_ang(const int dec = 30); //implementa pixerizacao
  void load_pix_file(const char *pixfile);
  void healpix_pixelization(const int dec); //pixelizacao do HEALpix
  void faz_super(void);
  void refaz_super(void);
  void faz_super(surface &Image);
  void faz_super2(surface &Image);
  void inline inc(double new_inc){par[0][0]=new_inc;};
  void inline M1(double nm1)
  {
    par[0][1]=nm1;
    roche_par.q=(par[0][2]/par[0][1]);
  };
  void inline M2(double nm2)
  {
    par[0][2]=nm2;
    roche_par.q=(par[0][2]/par[0][1]);
  };
  double inline inc(void){return par[0][0];};
  double inline M1(void){return par[0][1];};
  double inline M2(void){return par[0][2];};
  double inline K2(void){
	  using namespace Fis_const;
	  double ww = 1./(par[0][3]*3600.); //frequencia angular
	  //sep. orbital
	  double a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4.0*M_PI*M_PI*ww*ww), 1./3.);
	  double r2 = 1./(par[0][2]/par[0][1]+1.);//raio orbital da sec. em unidade de a
	  return 2.*M_PI*ww*a*r2/1e5*sin(par[0][0]);
	  
  };
	
  int inline dec(void){return __Dec__;};
  double inline ffac(void){
	return roche_par.f_fac;
  };
	void inline ffac(double new_ffac){
		roche_par.f_fac = new_ffac;
		refaz_super();
	};
	double R2_rsun(void);
	double inline oculta(vector<double>& p,vector<double> &opos){
		//Calcula se lobo de roche oculta pixel na posicao x,y,z no vetor pos
			
			// rotacionando ponto p em fase e em inclinacao
			double x = p[0]-1.0;
			double y = p[1];
			double z = p[2];
			double yp = y*opos[3]-x*opos[2];
			double xp = x*opos[3]+y*opos[2];
			double zp = z*opos[1]-xp*opos[0];
			
			if (xp >0 && sqrt(yp*yp + zp*zp) < ache_r(M_PI/2., 0.0)) {
				return 0;
			}
			
			x = p[3];
			y = p[1];
			z = p[2];
			yp = y*opos[3]-x*opos[2];
			xp = x*opos[3]+y*opos[2];
			zp = z*opos[1]-xp*opos[0];
			//cout << yp << " " << zp << " " << xp << endl;	
			
			
			return 1.0;
		
		
	};

};



//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+  


class SpotLight: public Light_Curve {
  static const int N_PAR_ORBITAL = 3; // numero de parametros orbitais
  // descricao:
  // 0 - i (radianos)
  // 1 - M_1
  // 2 - M_2
  static const int N_PAR_EMISS = 5; // parametros de emissao
  // descricao:
  // 0 - fluxo maximo
  // 1 - = "-1" (lei de obscurecimento de borda) 
  // 2 - direcao de maxima emissao (em graus)
  // 3 - Raio onde se localiza o spot (cood. esfericas)
  // 4 - theta (angulo em relacao ao eixo x, em graus)
  // 5 - dtheta (largura a meia altura de uma gaussiana para o spot, em graus)

  double Fi_Roche;
  double RL1;

  void lagrangeana(void);
  double Fi(double x, double y, double z); 
  void v_ang(const int dec = 30);
  void faz_super(void);
  double RaioLob1(double v0, double u0);

 public:
  SpotLight(const int dec, const vector<vector<double> >  &data_par);
  int N_Imax; // posicao no vetor do centro do spot
  double flux_ec(double fase);
  double spot_ec(vector<double> &, vector<double> &);

};

class Doppler: public Light_Curve {
  static const int N_PAR_ORBITAL = 4; // numero de parametros orbitais
  // descricao:
  // 0 - i
  // 1 - M_1
  // 2 - M_2
  // 3 - P_orb
  // 4 - V_max (km/s, V_min = -V_max)
  static const int N_PAR_ATMOSF = 0; // numero de parametros atmosfericos
  // descricao:
  // 0 - brilho da estrela em fase 0		       
  // 1 - coef. para obscurecimento gravitacional
  // 2 - numero de spots na superficie da estrela      
  // 3 - numero de coef. para obscurecimento de borda  
  // 4 - Temperatura no polo da estrela                

  struct par_list roche_par;
  double rl1; // raio do ponto lagrangiano interno
  void lagrangeana(void);
  void pixelization(const int SIZE); //pixerizacao...
  void faz_super(void);
  void faz_super(const surface &Image);

 public:
  Doppler(const int SIZE, const vector<vector<double> > &data_par);
  Doppler(const surface dmap, const vector<vector<double> > &data_par);
  //Doppler(char *pixfile, const vector<vector<double> > &data_par);
  //Doppler(const surface &StarSurf, const vector<vector<double> > &data_par);
  void Doppler_Image(surface &Image);
  //void GL_pix(const int dec, surface &Image);
  //double flux2(void) {return F20;};
 
};

//template <typename T> double Light_Curve::flux(const double fase, double (&ec_fnc)
//    (...),T ec_par)
template <typename T> double Light_Curve::flux(const double fase, T &ec_par, double (&ec_fnc)(vector<double>&,vector<double>&,T&))
{
	double n_s, n_s2, u1, lt_o, rlob, x1, y1, z1, z2, xp, yp, zp, f2, f1;
	double x, y, z, t;
	int k, nobs = int(par[1][2]);
	
	k = 0;
	lt_o = 0.;
	//u1 = fase*2.*M_PI;
	u1 = fase*2.*M_PI;
	double ffase = fase;
	z1 = cos(par[0][0]);
	z2 = sin(par[0][0]);
	f1 = sin(u1);
	f2 = cos(u1);
	
	x1 = z2*f2; //sin(par[0][0])*cos(u1);
	y1 = f1*z2; //sin(u1)*sin(par[0][0]);
	
	/*
	 x1 = sin(par[0][0])*cos(u1);
	 y1 = sin(u1)*sin(par[0][0]);
	 z1 = cos(par[0][0]);
	 z2 = sin(par[0][0]);
	 f1 = sin(u1);
	 f2 = cos(u1);
	 */
	//printf("%f", rl1-lagrangeana(1/d[1]));
	
	for (int i = 0; i < star.theta.size(); i++){
		
		n_s = star.nx[i]*x1+star.ny[i]*y1+star.nz[i]*z1;
		//    n_s2 = star.px[i]*x1+star.py[i]*y1+star.pz[i]*z1;
		
		n_s = n_s/sqrt(star.nx[i]*star.nx[i]+star.ny[i]*star.ny[i]+star.nz[i]*star.nz[i]);
		
		if(n_s > 0){ 
			// passo posicao (x,y,z) do pixel e z1,z2,f1,f2 para eclipse... ele que se vire com isso...
			vector<double> pos(4);
			pos[0] = star.px[i]; //-1./(par[0][2]/par[0][1]+1.);
			pos[1] = star.py[i];
			pos[2] = star.pz[i];
			pos[3] = star.px[i]-1./(par[0][2]/par[0][1]+1.);    
			vector<double> opos(4);
			opos[0] = z1;
			opos[1] = z2;
			opos[2] = f1;
			opos[3] = f2;
			
			lt_o += star.I[i]*n_s*obscurecimento_de_borda(n_s,nobs)*ec_fnc(pos,opos,ec_par);
		}
	}
	// Aplicando fator de escala
	
	if(par[1][0] > 0) lt_o *= par[1][0];  
	
	//  fazer isso em outro lugar!!!
	// lt_o = lt_o+d[4]; // d[4] = brilho do disco..
	
	return lt_o;
	
}

#endif